/*
 ============================================================================
 Name        : mBuild_inhomogenous.cpp
 Author      : Yi He
 Version     :
 Copyright   : A free distributable Academic program
 Description : The main program to build a membrane.
 ============================================================================
 */

#include <iostream>
#include "time.h"
#include "mbuild.h"

using namespace std;
void find_file(char* inFile,char* type,char* readfile,char* outfile);
void build_diel(char* fname,ToroidalPore &pore);
void build_charge(char* fname,ToroidalPore &pore);
void build_kappa(char* fname,ToroidalPore &pore);
void build_phi(char* fname,ToroidalPore &pore);
void build_shape(char* fname,ToroidalPore &pore);
void setup_pore(ToroidalPore&);

int main(void) {
	/* Input a file name to guess the other files*/

	  char inFile[256];
	  cout<<"Input a file name to get all other files:\n";
	  cin>>inFile;

	  //Setup pore
	  ToroidalPore pore;
	  setup_pore(pore);
	  build_diel(inFile,pore);
	  build_charge(inFile,pore);
	  build_kappa(inFile,pore);
	  build_phi(inFile,pore);
	  build_shape(inFile,pore);
	  return 0;
}

void setup_pore(ToroidalPore &pore)
{
  double r0,k0,xl,t,xh,yl,yh,zl,zh,cutoff;
  int ncyc;
  double r_up,r_bottom;
  cout<<"Set up Toroidal pore size(set R0=0 and K0=0 if don't want to have a toroidal pore).\n";
  cout<<"Pore radius R0=";
  cin>>r0;
  cout<<"K0=";
  cin>>k0;
  cout<<"Membrane thickness T=";
  cin>>t;
  cout<<"Set up hole(for cylindrical pore)\n";
  cout<<"R bottom=";
  cin>>r_bottom;
  cout<<"R up=";
  cin>>r_up;
  pore.setup(r0,k0,t);
  CylindricalHole hole=CylindricalHole(r_up,r_bottom,t+10*cutoff,-t/2-5*cutoff);
  if (r_up >0 && r_bottom >0) {
      pore.has_hole(&hole);
   }
  cout<<"Build region:\n";
  cout<<"X min=";
  cin>>xl;
  cout<<"X max=";
  cin>>xh;
  cout<<"Y min=";
  cin>>yl;
  cout<<"Y max=";
  cin>>yh;
  cout<<"Z min=";
  cin>>zl;
  cout<<"Z max=";
  cin>>zh;
  cout<<"Setup integration parameter:\n";
  cout<<"Integration cutoff:";
  cin>>cutoff;
  cout<<"Integration steps:";
  cin>>ncyc;
  pore.set_region(xl,xh,yl,yh,zl,zh);
  pore.set_cutoff(cutoff);
  pore.set_integration_cyc(ncyc);
}

void build_diel(char* fname,ToroidalPore &pore)
{
  double memb,head,water,width,diel_p;
  cout<<"Build dielectric maps.\n";
  //Setup Switch for dielectric constant;
  cout<<"Setup dielectric distribution:\n";
  cout<<"Dielectric constant in membrane:";
  cin>>memb;
  cout<<"Dielectric constant in Head group region:";
  cin>>head;
  cout<<"Dielectric constant in Water:";
  cin>>water;
  cout<<"Head group region  width:";
  cin>>width;
  cout<<"Solute dielectric constant(to skip):";
  cin>>diel_p;
  DielectricConstant c;
  c.set_diels(water,head,memb,diel_p);
  c.set_headgroup_width(width);
  pore.set_distribution(&c);
  GridSpace grid,reference;
  char input[256];
  char output[256];
  find_file(fname,"dielx",input,output);
  cout<<"Reading "<<input<<endl;
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);
  find_file(fname,"diely",input,output);
  cout<<"Reading "<<input<<endl;
  t1=clock();
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
  t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  find_file(fname,"dielz",input,output);
  cout<<"Reading "<<input<<endl;
  t1=clock();
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  pore.fill(&grid,&reference);
  cout<<"Creating hole\n";
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
  t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
}

void build_charge(char *fname,ToroidalPore &pore)
{
  double offset,area,anfr,width;
  cout<<"Setup charge distribution.\n";
  GridSpace grid,reference;
  int n=2;
  char read[256];
  char output[256];
  ChargeDensity c=ChargeDensity(n);
  QuadraticChargeDistribution d;
  double a,b;
  a=(20*pore.r0*pore.k0/pore.T+6*pore.k0*pore.k0/pore.T-5*pore.r0-3*pore.k0)/(10*pore.r0+2*pore.k0);
  b=4*(1-a)/(pore.T*pore.T);
  d.setup(a,b,pore.r0,pore.k0,pore.T);
  find_file(fname,"charge",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Head group has N layers of charge:";
  cin>>n;
  c.set_nlayer(n);
  for(int i=1; i<=n;i++)
  {
  cout<<"Layer  "<< i<<":\n";
  cout<<"Charge per lipid:";
  cin>>anfr;
  cout<<"Area per lipid:";
  cin>>area;
  cout<<"Offset of charged plane:";
  cin>>offset;
  cout<<"Gaussian distribution width:";
  cin>>width;
  c.set_layer(i-1,offset,width,anfr/area);
  }
  pore.set_distribution(&c,&d);
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
}

void build_kappa(char* fname,ToroidalPore &pore)
{
  double width;
  cout<<"Setup kappa\n";
  cout<<"Set ion inacessible region width:";
  cin>>width;
  Kappa k;
  k.set_headgroup_width(width);
  pore.set_distribution(&k);
  GridSpace grid,reference;
  char read[256];
  char output[256];
  find_file(fname,"kappa",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);

}

void build_phi(char* inFile,ToroidalPore &pore)
{
  Gouy g;
  double temp,offset,conc,z,area,anfr;
  cout<<"Generating phi based on gouy-chapman\n";
  cout<<"Temp=";
  cin>>temp;
  g.set_temp(temp);
  cout<<"Offset=";
  cin>>offset;
  g.set_offset(offset);
  cout<<"Ion concentration=";
  cin>>conc;
  cout<<"Valence=";
  cin>>z;
  g.set_ion(conc,z);
  cout<<"Area per Lipid=";
  cin>>area;
  cout<<"Fraction of anionic lipids=";
  cin>>anfr;
  g.set_lipid(area,anfr,1.0);
  pore.set_distribution(&g);
  GridSpace grid;
  char read[256];
  char output[256];
  find_file(inFile,"charge",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  cout<<"Finished.\n";
  grid.clear(0);
  find_file(inFile,"phi",read,output);
  cout<<"Building "<<output<<" ...\n";
  clock_t t1=clock();
  pore.fill(&grid,&grid);
  pore.create_hole(&grid,&grid);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);

}
void build_shape(char* fname,ToroidalPore &pore)
{
   GridSpace grid;
   char read[256];
   char output[256];
   IMM1Transition shape;
   double R0,K0,T;
   R0=pore.r0;
   K0=pore.k0;
   T=pore.T;
   shape.set_nsmth(10);
   shape.set_memb_thickness(T);
   pore.set_distribution(&shape);
   find_file(fname,"charge",read,output);
   DxFile::read_dx(read,&grid);
   grid.clear(0);
   find_file(fname,"shape",read,output);
   cout<<"Building "<<read <<endl;
   clock_t t1=clock();
   pore.fill(&grid,&grid);
   pore.create_hole(&grid,&grid);
   clock_t t2=clock();
   cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
   DxFile::write_dx(output,&grid);

}
void find_file(char* inFile,char* type,char* readfile,char* outfile)
{
  int l=strlen(inFile);
  char mid[256];
  strncpy(mid,&inFile[5],l-8);
  int l1=strlen(type);
  strncpy(readfile,type,l1);
  readfile[l1]='\0';
  strcat(readfile,mid);
  readfile[l+l1-8]='\0';
  strcpy(outfile,readfile);
  strcat(readfile,".dx");
  outfile[l+l1-8]='\0';
  strcat(outfile,"m.dx");
}
